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ABSTRACT 

■ We examine the validity of the hydrostatic equilibrium (HSE) assumption for galaxy clusters using 
one of the highest-resolution cosmological hydrodynamical simulations. We define and evaluate several 
effective mass terms corresponding to the Euler equations of the gas dynamics, and quantify the degree 

■ of the validity of HSE in terms of the mass estimate. We find that the mass estimated under the HSE 
assumption (the HSE mass) deviates from the true mass by up to ~ 30 %. This level of departure from 
HSE is consistent with the previous claims, but our physical interpretation is rather different. We 
demonstrate that the inertial term in the Euler equations makes a negligible contribution to the total 
mass, and the overall gravity of the cluster is balanced by the thermal gas pressure gradient and the 
gas acceleration term. Indeed the deviation from the HSE mass is well explained by the acceleration 
term at almost all radii. We also clarify the confusion of previous work due to the inappropriate 
application of the Jeans equations in considering the validity of HSE from the gas dynamics extracted 

, from cosmological hydrodynamical simulations. 

U Subject headings: cosmology: theory - galaxies: clusters: general - methods: numerical - X-rays: 
galaxies: clusters 

^ ■ 

a.: 

Q ■ 1. INTRODUCTION 

Clusters of galaxies are important sources of various cosmological and astrophysical infor mation especially o n the 
^2 \ formation history of the large scale structure and the estimates of cosmological parameters (jAllen et al.ll201lT for a 
recent review). Among others, mass of clusters is one of the most fundamental quantities in virtually all studies. 
The most conventional method is based on X-ray observations of the intracluster medium (ICM) combined with 
the assumption that the ICM is in hydrostatic equilibrium (HSE) with the total gravity of the cluster (We call the 
mass estimated by this method "the HSE mass"). It is unlikely, however, that the HSE assumption strictly holds, 
O^l ' especially for unrelaxed clusters given their on-going dynamical evolution. Therefore it is important to examine the 
C [ validity of the HSE assumption, which has been mostly assumed just for simplicity. The quantitative analysis of its 
~* . validity and limitation is directly related to the applicability to the future scientific opportunities, including upcoming 
-ray missions such as extended Rontgen Survey with Imaging Telescope Array 6 (eROSITA) and ASTRO-H 7 , and 
£SJ ' observations of the Sunyaev-Zel'dovich effect performed by Atacama Cosmology Telescope 8 (ACT) and South Pole 
O ■ Telescope 9 (SPT). 

The validity of the HSE assumption for observed clusters may be examined in a straightforward fashion by comparison 
of the HSE mass with the cluster mass estimated by other methods. In this respect, gravitational lensing is particularly 
suited because it directly probes the total gravitational mass without any assumption on the dynamical state of dark 
matter. On the other hand, the lensing observations require a high angular resolution of the background galaxy 
images and are feasible only for a limited number of clusters located at z < 0.5. In addition, the estimated lensing 
mass correspond to the cylindrical mass along the line of sight, and may include an ex tra contribution not associated 
with the cluster itself. Previous studies (e.g.. iMahdavi et al.ll2008t IZhang et al.ll2008L for recent ones) show that the 
HSE mass is smaller by approximately 20 percent on average than the lensing mass, suggesting either HSE or lensing, 
or even both, should be systematically biased. 

Another method to examine the validity of the HSE assumption that we pursue in this paper is to use numerical 
simulations, which enable us to make a detailed and critical comparison of the simulated data against the model 
prediction. Therefore we can locate the origin of systematic bias, if any, of the HSE assumption. This is useful because 
we may be able to apply the correction to the observational data eventually. 

Of course there are a number of previous studies of HSE using sim ulated clusters, but their results d o not seem to 
be converged. For instance, let us focus on a couple of recent papers (|Fang et al.ll2009t lLau et aLll2009D that studied 
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system atic errors in the HSE mass using the same set of 16 clusters simulated by iNagai et al.l (|2007l ). iFang et all 
(2009) analyzed the gas particle data on the basis of t he Euler equatio ns, and evaluated the effective mass terms 
corresponding to several different terms in the equations. lLau et ahl (|2009T ) performed basically the same analysis, but 
used the Jeans equations instead of the Euler equations despite the fact that they considered the gas particles in the 
simulated clusters. Both reached the similar conclusion that the HSE mass underestimates the true mass of clusters 
systematically by ~ 10-20 %. Nevertheless their physical interpretations of the origin of the bias are very different; 
IFang et ail ()2009[ ) claimed that th e coherent rotation of gas plays a significant role as an additional support against 
the gravity, while lLau et al.l (|2009l ) concluded that the random gas motion is responsible for the departure from HSE, 
and the gas rotation makes a relatively negligible contribution. 

The purpose of the present paper is to clarify the theoretical method to examine HSE of the simulated clusters, and 
then to revisit its validity using a high-resolution h ydrodynami c al sim ulati on by ICenl (120121) . In particular, we compare 



the two different analysis formulations adopted by IFang et all (|2009h and lLau et all ( 20091) , and argue tha t the Euler 
equat ions, rather than the Jeans equations modified by a gas pressure gradient term (jRasia et al.1 12004 lLau et ail 
2009), should be used in analyzing the gas dynamics. 

The rest of the paper is organized as follows. Section 2 formulates our method to evaluate the validity of the HSE 
assumption. We pay a particular emphasis on the comparison between the Euler and Jeans equations in analyzing 
the gas dynamics. The derivation of those sets of equations starting from the Boltzmann equations is summarized in 
Appendix A. Section 3 briefly describes the simulated cluster, and then presents our analysis results of the validity of 
HSE. Finally our conclusion is summarized in Section 4. The analysis based on the Jeans equations using dark matter 
particles is shown in Appendix B. 

2. THEORETICAL FORMULATION 
2.1. Method to Examine the Validity of HSE based on the Euler Equations 
Our analysis method to discuss the validity of HSE using gas in a simulated cluster is based on the Euler equations: 

dv I 

— + {v-V)v = Vp-V0, (1) 

at p gas 

where <f> is the gravitational potential, and p gas , v and p are the d ensity, veloci t y and pres sure of gas . Whil e the Jeans 
equations do not describe the gas dynamics in their orig inal form. iRasia et al.l f2004) and lLau et all (|2009fl added the 
gas pressure gradient term to the Jeans equations, and adopted the resulting equations in analyzing simulated clusters. 
We argue in the next subsection that this is not justified, and present the result based on the Jeans equations but 
using collisionless dark matter particles in Appendix B, for reference. 
We define the total mass M tot of a cluster inside a volume V as 

Mtot = / d 3 x ptot, (2) 
Jv 

where ptot is the total density of the cluster. For the simulated cluster considered throughout this paper, p to t consists 
of densities of gas, dark matter and stars, i.e., ptot = Pgas + Pdm + Astar- The total mass can be rewritten in terms of 
p and v using Poisson's equation and Gauss's theorem: 



Mtot = tt? / dS ■ Vcf> = / dS ■ 
4irG J dv AttG J dv 



1 w , ™ dv 
v>- [v -\/)v - — 

Pgas at 



(3) 



where dV is the surface surrounding the volume V and we have used equation (JTJ) in the second equality. Now the 
total mass is evaluated by the gas quantities alone, without any knowledge on dark matter and stars. This is why the 
present method is applicable, in principle, to the X-ray data of galaxy clusters. 
If we adopt a spherical surface as dV , the total mass can be decomposed into the following four effective mass terms: 

M tot = Mtherro + M rot + Mstrcam + M acce l, (4) 

where 

If 1 dp 

M therm = -— / dS— ' (5) 



4ttG 

lev Pgas dr 

M Iot = -±-[ dS V -^, (6) 
4ttG J 9v r 



jWst.rHa.Tn — : — / dS 

dV 



4ttG 

and 



dv r vg dv r Vu, dv r 



dr r 89 r sin 6 dip 



(7) 



4ttG J dV dt 



dS (8) 



We emphasize here that the above set of equations does not assume spherical symmetry of the system; we just take a 
spherical surface as the integral surface and write down equation Q in spherical coordinates. 



Validity of HSE 



3 



The first term, Mtherm, originates from the thermal pressure gradient of gas. If the gas motion is negligible, M t herm 
should be equal to the total mass M tot , and thus regarded as a cluster mass estimated under the HSE assumption. 
In other words, the difference between M to t and Mtherm is a quantitative measure of the departure from the HSE 
assumption. 

The incrtial term (v ■ V)v in the Euler equations reduces to M rot and M s t rC am- If there exists a coherent rotational 
motion around the center of the cluster, Af rot can be interpreted as the centrifugal force term. Without such a motion, 
however, the local tangential velocity of different directions at different locations on the sphere could make M rot 
significantly large. During the course of the cluster evolution, gas generally falls toward the center of the cluster with 
larger streaming speed with increasing radius from the center. In this case, M st[cllm becomes negative and cannot be 
neglected. On the other hand, it becomes positive and/or negligible, especially in the innermost region where the gas 
velocity is more randomized than that in the outer regions. 

Finally the acceleration term, M acco i, corresponds to the —dv r /dt term in the Euler equations, and becomes posi- 
tive/negative when gas is decelerating/accelerating. 

All the mass terms are invariant with respect to the choice of the axis of the spherical coordinates, but are not 
necessarily positive. Note also that M rot and M s t r eam, corresponding to the inertial term, are not invariant with 
respect to the Galilean transformation. Thus we evaluate those in the center-of-mass frame of the entire simulated 
cluster. 

2.2. Comparison with Analysis Methods Adopted by Previous Work 

The set of basic equations that we adopt in this paper is essentially identical to that of iFang et al.l (|2009t l. except 
the fact that they interpreted the difference between M to t and Mtherm + M rot + M s t r eam as "turbulent gas motion" 
while we call it the acceleration term M accc i as directly implied from equation ([1]), and evaluate it from the residual, 

M acC cl = Mtot — Mtherm ~ M rot — M s t r eam- 

We are not sure why they ascribed the term to the turbulent motion. It is true that part of the gas acceleration 
would be due to the turbulent gas motion, but not entirely. Furthermore the numerical simulation does not include 
any physical processes directly related to the turbulent motion. Even if the turbulent motion might be important 
for real clusters, it should come from some physics below the subgrid scales that cannot be properly resolved in the 
numerical simulation. Effects of gas random motion above the resolved scales should be included in M rot and Mstream - 
We note, however, that the different interpretation of M aC coi does not affect at all the conclusion of IFang et al.l (2009) 
that gas rotation is the most important term to describe the origin of departure from HSE in their simulated clusters. 
As we will show below, this is not consistent with our result. 

In previous literature, the Jeans equations are sometimes used in analyzing the gas motion in s imulated clusters . 
For that purpo se, a thermal pressure gradient term is added by hand to the basic equations (e.g.. lRasia et al.ll2004t 
lLau et aDl2009t ). Assuming the steady state, i.e., dv/dt = 0, the Jeans equations in r-direction (|A20|) is now replaced 
by the following equation: 



d ve d v v d 



dr r d6 r sin 9 dtp 



V r 



1 
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1 djp^afg) 1 d(p gas 

a rip) 
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, cot# 
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(9) 



r pga S dr r p gas dr dr ' 

where afj denotes the ij-component of the velocity dispersion tensor. 

For example, lLau et all J2009) converted each term of equation ([9]) to the corresponding mass term; see their 
equations (6) to (11). Then they interpreted that the first two terms on the right hand side of the above equation 
originate from "random gas motion" . We suspect that the implicit assumption underlying equation ^ is that the 
ICM consists of two distinct gas components; the thermal gas and the unthermalized one such as the cold gas accreting 
along with galaxies into clusters. While it may be reasonable to take into account both components separately when 
dealing with real galaxy clusters, it remains to be justified if the extra pressure gradient term can be added by hand to 
the conventional Jeans equations assuming that the two components share the same density p gas (In Appendix A we 
show that diagonal components of the velocity dispersion tensor in the Jeans equations correspond to thermal pressure 
in the Euler equations). It should also be recalled that the current numerical simulations are performed entirely on 
the basis of the Euler equations. Therefore, for the sake of consistency, it is appropriate to use the Euler equations, 
instead of equation (J9j) , in evaluating the mass terms of the simulated galaxy clusters. The proper treatment of the 
unthermalized gas in numerical simulations is certainly an important issue, but is beyond the scope of the present 
paper. 

3. APPLICATION TO A SIMULATED CLUSTER 

We apply the analysis method in the previous sect i on and quantitatively examine the validity of the HSE assumption. 
For that purpose, w e use a simulated cluster oflCen] (l2012t) . It is simulated with an Eulerian ad aptive mesh refinement 
(AMR) code, Enzo (|Brvanlll999HBrvan fc Normarifl999T-TO'Shea et "aTll2004t [Joung et al.ll2009t ). Although our current 
analysis exploits three-dimensional profiles of the cluster which cannot be observed in reality, we are interested here 
in the validity of HSE itself. The application to the X-ray observation will be discussed separately elsewhere. 
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3.1. Hydrodynamical Simulation 

We briefly describe the main features of the current simulation, and refer readers to ICenl ()2012f ) for more detail. 
The simulation was run first with a low resolution mode in a periodic box of 120 hr x Mpc on a side. Then a region 
centered on a cluster with a mass of ~ 2 x 10 14 /i _1 M Q was resimulated with a higher resolution in an adaptively 
refined manner. The size of the refined region is 21 x 24 x 20(/i~ 1 Mpc) 3 . The mean interparticle separation and the 
dark matter particle mass in the refined region are 117 h~ 1 kpc comoving and 1.07 x 10 8 h~ 1 M ( ? ) , respectively. 

Star particles are created according to the prescriptio n of ICen fc Ostriker (1992). Their typical mass is ~ 1O 6 M . 
The simul ation includes a m etagalactic UV background (|Haardt fc Madaulll996F). shieldin g of UV radiation by neutral 
hydrogen ()Cen et al I2005D and meta llicity-dependent radiative cooling (|Cen et al"1 ll995). While supernova feedback 
is modeled following ICen et all ()2005[ ). AGN feedback is not included in this simulation. The cosmological parameters 
used in this simulation ar e (fl ft , » m , Q A , fe, n, g R ) = (0.046, 0.28, 0.72, 0.70, 0.96, 0.82), following the WMAP7- 
normalized ACDM model (|Komatsu et al.ll201lT) . 

Then a cluster is identified and the cubic box of a side of 3.8/i _1 Mpc surrounding the entire cluster is extracted 
from the simulation data. The dark matter and stars are represented by particles, and the temperature and density 
of gas are given on the 520 3 grids (the grid length is 7.324ft~ 1 kpc). 

The radius r$oo of the cluster is ~ 640ft. -1 kpc (r$oo is defined so that the mean density inside rsoo is 500 times the 
critical density of the universe) . The center-of-mass velocity of the cluster within rsoo is set to vanish. The total mass 
-W500 within r5oo is ~ 2 x 10 M Q . The average ICM temperature at r5oo is ~ 2 keV, and the circular speed there is 

W500 = y/GM 5m /r 5 oo ~ 1000 km s _1 . 

Projected surface densities of gas, dark matter and stars on x-z plane are plotted in the left panels of Figure [1] The 
right panels of Figured] show the three-dimensional view of the three components. The left (right) plots are color-coded 
according to the surface (space) densities normalized by the fraction of each component averaged over the box, fi^ 
(fc = gas, dark matter and stars). Note that the fraction Cl^ is different from the density parameter flk because the 
box is selected preferentially around the cluster. 

As is clear from FigurelU the gas distribution is smoother but very well traces the underlying dark matter distribution. 
In contract, stars are more significantly concentrated in high density regions, and exhibit numerous small clumps, most 
of which are not identified/resolved in the gas distribution. 

Figure [2] plots the radial density and mass profiles of the cluster. The stellar fraction in the inner region (r < 
200ft~ 1 kpc) is significantly higher than the typical observed value. This is a well-known common problem among 
current high-resolution cosmological simulations, and implies that some important baryon physics including high- 
energy phenomena and star formation is still missing in the simulation. We perform the analysis of cluster gas, 
assuming that this excessive star densities in the inner region does not affect our conclusions at outer radius. 

Figure [3] represents velocity fields in x-y, y-z and z-x planes passing through the center of the cluster. The red/blue 
arrows have negative/positive radial velocity, showing that the gas in the outer regions (r > Ih^ 1 Mpc) falls toward 
the center while its direction is randomized in the inner region. 

3.2. Results 

We evaluate the effective mass terms defined in Section 2 for the simulated cluster, which is plotted in Figure [SJ the 
left panel shows the mass profiles, while the right panel indicates their fractional contribution to the total mass within 
the radius. 

The total mass M to t(r) is computed by directly summing up all the dark matter and star particles and gas of grids 
within the sphere of r. The other terms, Mtherm> M TO % and M s t r eamj require the pressure and velocity fields evaluated 
at r. For that purpose, we use the density, velocity and temperature of gas defined at 520 3 original grid points, and 
first bin the cluster into 50 logarithmically equal radial intervals between 90-1900 h^ 1 kpc, 10 x 10 linearly equal 
angular intervals. 

Integrands of equations (0, ([6]) and ([7|) are calculated in each bin. Derivatives of the physical variables such as dp/dr 
are calculated as follows; first d/dx, d/dy and d/dz are computed from the difference of the adjacent the original 
Cartesian grid points. Then d/dr, d/dO and d/dip in our spherical coordinates are calculated applying the chain rule. 

In this way, Mtherm(^)) M rot (r) and M stT ea,m(r) are computed by integrating the corresponding integrands evaluated 
above. Finally we estimate M acco i simply from the residual of M acce i = M to t — M t h e rm — M TO t — M s t rcam , since we have 
the cluster data at z = alone. This estimate for M accG i may be different from the original defintion, i.e., equation (jSJ). 
Indeed when we attempted to compute M accc i directly from the box mentioned in Section 3, it turned out to be too 
small to obtain a correct gravitational potential for the entire cluster. Thus we go back to a larger simulation box of a 
side of 22.5 h^ 1 Mpc and the grid length of 29.34 h^ 1 kpc that encloses our cluster. Then we compute the gravitational 
potential using FFT to obtain the gas acceleration at each grid point. This enables us to directly calculate M acce i. 
Figure [4] is a comparison of M acce i's calculated by two methods. The directly calculated M acco i (magenta line) is in 
good agreement with M to t — Mtham — M Tot — M strC am (black line), although there is a large difference between the 
two within 200h~ 1 kpc. Also, we make sure that the sum Mtherm + Mrot + -^stream + -Maccei reproduces M to t within 
~ 2 % except for the innermost region (r < 200ft~ 1 kpc), where it deviates from M to t by up to ~ 9 %. Thus the 
estimation of Maccei by Af to t — M^ elm — M rot — M str cam is sufficiently good given the quoted errors of our conclusion 
below. Although it seems better to use M acC oi directly calculated from equation ((HJ, the grid size of the larger box is so 
coarse that we cannot take advantage of the high resolution of the simulations in this study. Therefore, we decided to 
use the smaller box explained in Section 3 and M acC ci is calculated by M to t — -Mtherm — Mmt — -^stream i n the following 
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analysis. 

The left panel of Figure [5] implies that Mthcrm agrees with Mtot reasonably well. Each of the other three terms 
contributes less than 10 % of the total mass (the dotted curves correspond to the case in which each term becomes 
negative and its absolute value is plotted instead). 

In order to consider the validity of HSE more quantitatively, we plot the fractional contribution of each mass term in 
the right panel of Figure [5l The rotation term, M rot , is always positive (by definition) and contributes approximately 
10 % almost independently of radius. In contrast, the streaming velocity term, Mgtreamj is mostly negative, and varies 
a lot at different radial bins. As a result, the difference of the total mass M tot and the HSE mass M t herm is mostly 
explained by the acceleration term M acce i alone; compare the black and magenta curves in the right panel of Figure 
At r = rsoo and T2oo, the deviation from HSE in terms of the mass difference (Mthcrm — M to t) / M to t is about 10 %. 
Nevertheless the value significantly varies at different radii and it is safe to conclude that (Mthcrm — Mt t)/M to t ranges 
approximately 10-20 % at r < r2oo- Also there is no systematic trend of the validity of the HSE assumption as a function 
of radius. Even though the reliability of the simulation is suspicious for r < 200/i _1 kpc due to the excessive stellar 
concentration (Section 3.1), (Mtherm — M t ot)/M to t fluctuates between —10 % and +25 % for 300ft- _1 kpc < r < r 5 oo- 
Thus there is no guarantee that HSE becomes a better approximation toward the inner central region. 

Physically speaking, M rot + M stre am + M acce i(= M to t — M t h C rm) corresponds to the term intergrating the Lagrangian 
derivative of the gas velocity over the sphere. Therefore the fact that it is small compared with Mtherm and M to t is 
simply translated into the condition of HSE that the gas acceleration from a Lagrangian point of view is negligible 
compared with the pressure gradient and the total gravity. 

Since we analyze a single simulated cluster, it is not clear to what extent our interpretation that M an cci determines 
the departure from HSE holds in general. Thus we use a different set of SPH simulation clusters (|Dolag et al.ll2009l) 
kindly provided by Klaus Dolag. For these clusters, we assume spherical symmetry and do not divide the spherical 
surface at a given radius just for simplicity. Klaus Dolag also provides us with M accc i directly computed from the 
acceleration data. We find that in most cases |M acce i| is larger than M rot and |M s t r0 am| especially where M t herm 
deviates from M to t. We also confirm that Mtherm + M rot + M s t rea m + M acce i reproduces M to t within ~ 5 % for most 
regions (See Appendix C). 

The above results basically support our conclusion for the cluster from the AMR simulations, but we do not have 
a statistical discussion including clusters from the SPH simulations because of the differences in resolutions and 
simulation methods. Instead, we divide the AMR cluster into two regions; upper and lower hemispheres with respect 
to the x-y plane. Then we duplicate each hemisphere into one cluster. We call the synthetic cluster constructed from 
the z > (z < 0) hemisphere "z+" ( "z-" ) . Although these clusters are of course not independent of the original cluster 
and we cannot make a statistical argument on their properties, we can briefly look at the effects of substructures or 
inhomogeneity of temperature and velocity field. 

We repeat the same analysis on these two synthetic clusters, and the results are plotted in Figure El The amplitudes 
of the different terms vary significantly between the two clusters, and the degree of the validity of HSE is also very 
different. Nevertheless the generic trend is clear; the relation of M accc i w M to t — Mtherm holds almost independently 
of r. 

It is not clear, however, why the two hemispheres have so different values of (Af to t — Mtherm)/Mt t; HSE holds very 
well for "z+" , while it is not the case for "z-" . The visual inspection of Figure[T]does not reveal any significant difference 
between the two. It may be because some local concentrations of dark matter enhance the acceleration/deceleration 
of gas, and influence the overall non-sphericity of the gas density. Thus the analysis taking account of the ellipticity 
may provide a deeper insight on the validity of HSE, but is beyond the scope of the present paper. 

4. CONCLUSION 

We have examined the validity of HSE that has been conventionally assumed in estimating the mass of galaxy 
clusters from X-ray observations. We use a simulated cluster and evaluate several mass terms directly corresponding 
to the Euler equations that govern the gas dynamics. We find that the mass estimated under the HSE assumption, 
Mtherm in the present study, deviates from the true mass M to t on average by ~ (10 — 20) % fractionally for r < r 2 oo- 
There is no clear tendency that the HSE becomes a better approximation toward the inner region. More importantly, 
we find that M to t — Mthcrm is nearly identical to M accc i , in other words, the validity of HSE is controlled by the amount 
of gas acceleration. This trend is confirmed by the separate analysis of the different hemispheres of the same simulated 
cluster. 

Our current analysis is limited to a single simulated cluster, but the ove rall conclusion tha t the HSE mass agree s 
with the total mass within (10-20)% is consistent with previous results bv lFang"eT al. (2009) and Lau et al.l {2009). 
Nevertheless the interpretation of the origin of the departure from HSE is very different. 
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(2009) concluded 



that the gas rotation term M rot makes a significant contribution and that Mthcrm + M rot well reproduce s the total 
mass, especially for relaxed clusters. It is not the case, however, for our simulated cluster at least. Similarly lLau et al.l 
(2009) found the similar degree of the departure from HSE, but they ascribed the discrepancy to the random gas 
motion. Their analysis, however, is based on the modification of the Jeans equations, which does not appear to be 
justified for the analysis of the gas dynamics, and thus their conclusion should be interpreted with caution. 

A relatively small systematic error of the HSE mass inferred from current numerical simulations may be partly 
ascribed to the assumptions inherent in the Euler equations, i.e., local thermal equilibrium and negligible viscosity 
(Appendix A) . This is supported by the fact that the error in the mass estimated from the random motion of collisionless 
particles tends to be much greater at large radii (Appendix B) because the relaxation time scale for collisionless particles 
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are appreciably longer than that for collisional gas. We should also note that the HSE mass can be influenced by other 
physical precesses that are not included in the num erical simulations, s uch as pressure support from micro-turbulence, 
the magnetic field, and accelerated particles (e.g.. lLagana et aLll2010fl . The neglected components mentioned above 
are closely linked with one another (e.g., viscosity can play a role in generating turbulence and the magnetic field 
can affect both thermalization and acceleration of gas particles) and will be investigated in the near future by X-ray 
missions such as NuSTAR 10 and ASTRO-H as well as by radio telescopes including EVLA 11 and LOFAR 12 . 

The present analysis is fairly idealized in a sense that the evaluation of all the mass terms has been done from the 
full three dimensional data of simulated clusters. In reality observational data of X-ray clusters are basically projected 
along the line of sight. Therefore additional uncertainties and systematic errors may become important as well. Also 
it is necessary to carry out the current analysis for a number of clusters in order to obtain the statistically robust 
conclusion. These issues will be discussed elsewhere. 
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APPENDIX 

RELATION BETWEEN THE EULER EQUATIONS AND THE JEANS EQUATIONS 

From a microscopic point of view, both the Euler equations and the Jeans equations can be derived from the 
Boltzmann equation under different assumptions. In the following, we explicitly compare the two sets of equations in 
both Cartesian and spherical coordinates. 



Cartesian Coordinates 

We define the distribution function / such that f(x,v,t)d 3 xd 3 v is the probability that a randomly chosen particle 
in the system lies in the phase space volume d 3 xd 3 v at position (x, v) and time t. The motion of such particles under 
the gravitational potential <fi is described by the Boltzmann equation: 

coll 



dt dx l dxi dv z \ 5t 



where the collision term on the right hand side takes account of collisions between particles. Note that v l (i = 1, 2, 3) 
represents a coordinate in the phase space and should not be confused with the velocity field at the spatial point x % . 
For simplicity, we assume that all particles have the same mass m in the following. 
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First, we consider a collisionless case with (6f/St) C 6i\ = 0. Multiplying equation (|A1I) by m and integrating it over 
the velocity space yield the continuity equation: 

where 

p(x,t) = J d 3 v mf(x,v,t) (A3) 
and we introduce the average over the velocity space: 

q(x,t) — - - y <i 3 u mq(x,v,t)f(x,v,i) (A4) 

for an arbitrary variable g such as v l . Multiplying equation (|A1[) by m!) J and integrating over the velocity space give 
the momentum equations: 

at + dxi ~ p dx^ (A5) 

where 

t] 3 = pv % vi = pa 2 ' ij + pv { v j (A6) 

and <7 2 ' u is the velocity dispersion tensor. Equations (| A2|) and (|A5I) reduce to the Jeans equations: 

dv l -.dtf ldtpa 2 ™) 86 

VvJ = — — (A7) 

at dxi p dxi dx t K ' 

Next, we consider a collisional case. Rigorous handling of the collisional term is rather complicated and simplified 
models are often used. A conventional one is the Bhartnagar-Gross-Krock (BGK) equation, which employs a linearized 
collisional term: 

+ _«£V = _/z£, (A8) 
dt dx l dxi dv l t ' 

where r is the relaxation time of the system considered, /o is the Maxwellian distribution function characterized by 
the local temperature T(x,t): 

-i 3/2 



fo(x,v,t) 



27rk B T(x,t) 



exp 



m(v — v(x, t)) 



- 1 



2k B T(x,t) 



(A9) 



and ks is the Boltzmann constant. If we assume that mean values of conservatives such as mass and momentum are 
the same as those in local thermal equilibrium, the collision term vanishes in the continuity and momentum equations. 
In local thermal equilibrium, pressure is defined from the diagonal components of afj by pSij = pa\^ and the dispersion 
tensor can be written as 

rf = P 5 ij +pv i v i . (A10) 



If we replace tj in equation (|A5[) with te and combine them with equation (|A2[) . we obtain the Euler equations: 

dv z ■ dv % 1 dp 36 , , x 

ot ox 3 p oxi axi 

The difference between the Euler and the Jeans equations resides only in the form of the dispersion tensor. 

If we retain the off-diagonal components of t he dispersion tensor, they can be interprete d as viscosity, and the 
equations reduce to the Navier-Stokes equations ()Choudhurilll998t IChapman fc Cowlindll97Cil) . 

Spherical coordinates 

One can rewrite the continuity and momentum equations in the previous section in general coordinates: 

^+V i 0w < ) = (A12) 



and 

' >U "" ] 1 V,-t« = -pg^V 3 6, (A13) 



dt 

by using the covariant derivative operator Vj . For spherical coordinates (x 1 = r, x 1 = 9, x 3 = cp), non-zero components 
of the metric tensor g lJ are 

011 = 1, 522 = r 2 , 533 = r 2 sin 2 9 (A14) 
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and the corresponding non-zero connection coefficients are 



r 1 

1 22 



p2 
1 33 



r 1 

1 33 



-r sin 



r 3 

1 13 



1 12 - 1 21 - "J 



p3 
1 31 



r 3 
1 23 



r^ 2 = COt ( 



The velocity vector is now given by v l = (v r ,vg/r, v v /r sm.9) 
In spherical coordinates, the continuity equation reads 

dp 1 d(j 2 pv. 
dt r' 



Setting t % - 



pv l v 3 



■ 2 dr r sin 9 89 

pg lJ gives the Euler equations: 



1 d(sm9pve) 1 d(pv v ) 
r sin 9 dip 



d_ 

dt 



8 



d vg d \ 
dr r 89 r sin 9 dip 



d Vg d 



a 



r 89 r sin 9 dtp 



dt r dr 

d d vg d 
dt 1 dr r 89 r sin 9 dip 



V r Vg — Vl COt ( 



v v d 



V r V v + VgV v COt 9 



1 dp d<j> 
p dr dr 

1 dp 1 dip 
pr d9 r d9 

1 dp 1 



pr sin 9 dp r sin 9 dip 



On the other hand, putting r 



ij - - 



pv l v 3 leads to the Jeans equations 



~d _ d 

dt r dr 






d ' 


V r 4 


1 


d{pa 2 r ) 


15(pa r 2 e ) 


1 




r 89 


r sin 6 


)8ip 
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dr 


r 96» 


r sin 6 


> dip 


+ \ - 


a lg ~ a l v 






+ a lb 


C01 


■>--£ 








~d _ d 

dt r dr 


+ ^- + 


v v 


d ' 


vg -f 


1 


d(p(T 2 g) 
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d(pv 2 g v )~ 


r 89 


r sin 6 


Wip_ 


p 


dr 


r d9 


r sin 6 


> dip 


1 / 9 

+ - 3a r 2 , - 
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o* v cot 9 - 


h v r vg 


- v 2 cot 6 


+ a 2 gg cot 9) = 


ld<t> 
~r~89 






' d d 

dt dr 


+ ^- + 
r 89 


v v 


d ' 


v v + 


1 




1 5(pa e 2 y ) 
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\ — 




r sm9 


dp _ 
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dr 


r <96» 


r sin ( 


3 dip 


f I (3<, + 


v r v v + vgi 


v cot t 


? + 2 


a 2 eip cot 9 


)= 1 
r sin 6 


dtp' 







SYSTEMATIC ERRORS IN MASS ESTIMATES FOR COLLISIONLESS SYSTEMS 
In a similar fashion to Section 2.1, we can compute the gravitational mass using the Jeans equations: 

dv i 

— + (17 • V)v = V(p dm <T 2 ) - V</>, 

dt pdm 



Aft, 



1 



4ttG 



1 , „ 9v 

V(p d m<T 2 ) - (17 • V)t7 - — 



(A15) 



(A16) 

(A17) 
(A18) 
(A19) 



(A20) 



(A21) 



(A22) 



(Bl) 
(B2) 



where v and er 2 are the velocity field of particles and the velocity dispersion tensor, respectively. We here represent 
the collisionless component by dark matter, but the same formulation is readily applicable to galaxies. We decompose 
the right hand side of equation (|B2I) into the following terms by means of equation (|A20[) : 



Aft, 



rand 



AL 



+ M TOt + M st 

ream ~\~ Af cross -\- M acce \. 



A/rand — 
Afaniso 

Afmt 



1 



4ttG 
1 

4^G 



dS 
av Pd 



1 d(p dm a 2 r ) 



dr 



dS 



2ol 



<)V 



AttG 



dS 



(B3) 
(B4) 

(B5) 

(B6) 



av 
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4ttG 



dS 



m ■ 



dr 



4nG 



dS 



1 d(p dm af, g ) 



Vg dv r 

~r~d9 
1 



dv r 



r sin dtp 
' 2 ) 

) 



cot( 



(9(9 



Pdm^ sin 9 dip 



1 /" ,c 



(B7) 
(B8) 
(B9) 



Physical interpretation of each mass term is as follows. The term M ran( j comes from the gradient of velocity dispersion 
in the r-direction and corresponds to M t horm for a collisional gas. The meaning of M rot , Af s t ream and Af accc i are similar 
to the corresponding terms for the collisional gas (equations (6) to (8)). The terms that have no counterpart in the 
Euler equations are M an ; so and M cross ; the former represents anisotropy of the velocity dispersion whereas the latter 
arises from the off-diagonal components of the velocity dispersion tenser and vanishes if velocities of different directions 
are uncorrelated. 

We apply the above formulation to dark matter particles in the simulated cluster described in Section 3 to quantify 
intrinsic systematic errors of the mass estimation using collisionless particles. Note that one can apply the same 
method to galaxies but with much larger impact of statistical errors. Therefore, we do not do so here because we are 
interested in intrinsic systematic errors independent of observational complexities. Each term is computed in a similar 
manner to the case of collisional gas described in Section 3. 

Figure [7] shows that the difference between M to t and M ran( j increases toward the outer envelope mainly owing to the 
presence of M an j S0 . This is because the relaxation timescale of collisionless particles is much longer than that of the 
collisional gas. Once this term is subtracted, Af to t — Af ran d — Af an i so closely matches M acco i whose absolute value is 
limited to within ~ 0.3M tot . The amount of Ajf acce i is similar to that for the collisional gas (Fig. [5]). The other mass 
terms such as Af cross are less important. 

The above results imply that proper account of velocity anisotropics is essential for the mass reconstruction using a 
collisionless component. We stress that Af an j so is irrelevant to the collisional fluid as long as local thermal equilibrium 
is established (Appendix A). 

VALIDITY OF HYDROSTATIC EQUILIBRIUM IN CLUSTERS FROM SMOOTHED PARTICLE HYDRODYNAMICS 

SIMULATION 

We have shown the behavior of the mass terms for the cluster from the AMR simulation in Section 4. Especially, the 
HSE mass deviates from the total mass by up to ~ 30 % and the difference is explained chiefly by the gas acceleration. 
In order to make sure if such behavior is common to simulated clusters, we analyze the other simulated clusters. 
Since our AMR simulation includes only a si ngle cluster, we an alyze clusters extracted from a smoothed particle 
hydrodyna mics (SPH) simulat ion performed bv lDolag et al.l (|2009f ) using GADGET-2 to confirm the results in Section 
4. We refer iDolag et al.l (|2009l ) for the detail of the simulation. 

We use five regions named gl, g72, gl542, g3344 and g914 extracted from the simulation. Each region has a cluster 
near the center and the cluster is labeled "a" (gla, g72a, gl542a, g3344a, g914a). For each cluster, we make 256 3 mesh 
data centered on the cluster's center-of-mass within the radius r 50 o. The size of the data box is determined so that it 
contains the entire sphere of radius 3r 50 o centered on the cluster's center-of-mass. 

We calculate the mass terms defined in Section 2 for the five simulated clusters. Since we have the gas acceleration 
data, we calculate Af acce i directly from the data without using Aftot — Afthcrm — Af ro t — Af s t r0 am. 

Figure [8] shows the result for the cluster g914a. This figure basically supports the results in Section 4 in that A/therm 
deviates from Af to t by ~ 30% at most and that Af acco i becomes large where the difference between Af to t and Af t hcrm is 



large. The fact that the sum M t hcrm + A/ r , 



Af accc i is approximately Af to t means the estimation of Af a 



by Aftot — Afthcrm — Af rot — Af s t rGam is good, as confirmed in Section 4 for the AMR cluster. 
Although not graphically shown, the mass terms for the other simulated clusters also exhibit similar behavior. 




1000 





-1000 1000 
x [ft. -1 kpc] 

Fig. 1. — Densities of gas (top), dark matter (middle) and stars (bottom) are shown in two ways for each. Left: Projected surface 
densities on x-z plane normalized by the fraction of each component in the box: logj^J dl p^Q^ 1 /(g cm -3 h~ 1 kpc)], where k is one of 
gas, dark matter and stars and Q k is the fraction of k component in the box. Right: Equal-density surfaces for log 10 [pj.f2^" 1 /(g cm -3 h —1 
kpc)]=— 28.0 (blue), —27.0 (green) and —26.0 (red). Densities are normalized in the same way as the left panel. 
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Fig. 2. — Radial profiles of densities (left) and masses (right) are shown for gas (red), dark matter (green) and stars (blue). The black 
line in the right panel shows the total gravitational mass; Mtot = Af gas + M dm + M s t ar • The analysis is performed on the 50 logarithmically 
equal radial bins. 
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Fig. 3. — Velocity fields in x-y (left), y-z (middle) and z-x (right) planes passing through the center of the cluster. The arrow is red if 
v r < 0, and blue if v r > 0. The length of the arrow is proportional to the magnitude of the velocity. An arrow with a speed of 1000 km 
s _1 is shown for reference. 



12 



Suto et al. 



10 



10 1 



15 



^tot ^therm ^rot ^stream 
directly calculated 





100 200 500 1000 

r [h _1 kpc] 



Fig. 4. — As comparison of the acceleration mass M acce i calculated directly from the acceleration data (magenta) and Mtot — Mtherm — 
M ro t — M s tream (black). Dotted line means that its sign is inverted. The analysis is performed on the 25 logarithmically equal radial bins. 




Fig. 5. — The effective mass terms in Equation J4} for the gas in the simulated cluster are shown in the left panel: Mtot (black), Mtherm 
(red), M r ot (green), M st ream (blue) and A/ accc i (magenta). Here iVf acce i is calculated by iVf acce i — Mtot — A^therm — iWrot — A/stream- Dotted 
line means that its sign is inverted. Ratios of mass terms to Mtot are shown in the right panel. The black line shows (Mtot — Mtherm) /Mtot 
and colored lines represent the same things as the left panel. The analysis is performed on the 50 logarithmically equal radial bins and 10 
linearly equal bins both in polar and azimuthal angles. 




Fig. 6. — Same as Figure l5l but for the clusters "z+" (top) and "z-" (bottom). 




Fig. 7. — The effective mass terms in Equation (IB3 I t for the dark matter in the simulated cluster are shown in the left panel: Mtot 
(black), M raIld (red), M rot (green), M st ream (blue) , M amso (cyan) , M croBB (orange) and M acce i (magenta) . Here M accc i is calculated by 
M acce i = Mtot — M ran( j — M rot — M st rcam — M amso — M croas . Dotted line means that its sign is inverted. Ratios of mass terms to Mtot are 
shown in the right panel; The black line shows (Mtot — M raru j)/Mtot and other mass terms are colored in the same colors as the left panel. 




Fig. 8. — The effective mass terms in Equation Q for the cluster g914a from the SPH simulation are shown in the left panel: Mtot 
(black), M thcrm (red), M rot (green), M at ream (blue), A^ amso (cyan), A<f cross (orange) and A/ accc ^ (magenta). It also shows the sum A4therm~i~ 
Mrot + Matream + M acce i in orange. Here M acce i is calculated using the acceleration data, not Af acco i = Aftot — M t i lerm — M ro t — M s t r eam- 
Ratios of mass terms to Mtot are shown in the right panel. The color-coding is the same as the left panel. The analysis is performed on 
the 60 logarithmically equal radial bins. 



